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The  purpose  of  this  research  has  been  to  further  develop  a  simple  efficient  grid 
generation  procedure  for  external  aerodynamics  aplicatiorft.  The  grid  generation 
scheme  is  based  on  solving  hyperbolic  partial  differential  equation  constraints  of  grid 
angularity  and  mesh  incremental  volumes.  The  grid  generation  scheme  has  been 
previously  used  in  two  dimensional  applications  to  generate  grids  about  smooth 
body  shapes.  The  main  thrust  of  this  AFOSR  supported  research  has  been  to 
extend  the  hyperbolic  partial  differential  equation  procedure  to  three  dimensional 
applications  and  to  study  ways  of  applying  the  procedure  to  body  shapes  that  have 
discontinuous  derivatives. 

The  main  part  of  this  report,  Part  I,  is  devoted  to  describing  the  three  dimensional 
hyperbolic  grid  generator.  This  Section  first  reviews  the  hyperbolic  grid  generation 
procedure  in  two  dimensions  and  then  describes  the  extension  to  three  dimensions.^ 
The  numerical  solution  algorithm,  mesh  control  functions  and  treatment  of  axis 
singularities  are  then  described.  Computed  meshes  and  even  are  calculated  flow 
fled  result  are  shown  to  help  evaluate  the  grids.  Part  I  of  this  report  is  essentially 
a  self-contained  technical  report  that  is  being  readied  for  submission  to  a  Journal. 
4-Part  II  of  this  report  is  both  brief  and  sketchy  in  its  presentation.  In  this  section 
we  describe  some  of  our  success  in  treating  bodies  with  sharp  edges  and  bodies 
that  are  exceptionally  concave. \The  hyperbolic  grid  generation  tends  to  propogate 
discontinuities  and  some  bodieslare  not  compatible  with  the  imposed  orthogonality 
and  volume  constraints  that  art  user  imposed.  When  this  happens,  the  hyperbolic 
grid  generator  breaks  dow^/  What  must  be  done  in  this  case  is  to  relax  these 
constraints.  We JiavediSdsome  success  doing  this,  but  this  process  remains  still  an 
art  form: 

^The  last  part  of  this  report  describes  a  flow  field  algorithm  development.  During 
the  course  of  this  research  we  had  some  considerable  interaction  with  AFWAL, 
and  at  one  point  became  ’side-tracked’  into  a  successful  approach  of  improving  the 
efficiency  of  our  general  implicit  Euler  aad  Navier-Stokes  code. i/P art  HI  of  this 
paper  contains  a  preliminary  paper  describing  this  development. 
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FOREWORD 


The  purpose  of  this  research  has  been  to  further  develop  a  simple  efficient  grid 
generation  procedure  for  external  aerodynamics  aplications.  The  grid  generation 
scheme  is  based  on  solving  hyperbolic  partial  differential  equation  constraints  of  grid 
angularity  and  mesh  incremental  volumes.  The  grid  generation  scheme  has  been 
previously  used  in  two  dimensional  applications  to  generate  grids  about  smooth 
body  shapes.  The  main  thrust  of  this  AFOSR  supported  research  has  been  to 
extend  the  hyperbolic  partial  differential  equation  procedure  to  three  dimensional 
applications  and  to  study  ways  of  applying  the  procedure  to  body  shapes  that  have 
discontinuous  derivatives. 

The  main  part  of  this  report,  Part  I,  is  devoted  to  describing  the  three  dimensional 
hyperbolic  grid  generator.  This  Section  first  reviews  the  hyperbolic  grid  generation 
procedure  in  two  dimensions  and  then  describes  the  extension  to  three  dimensions. 
The  numerical  solution  algorithm,  mesh  control  functions  and  treatment  of  axis 
singularities  are  then  described.  Computed  meshes  and  even  are  calculated  flow 
fied  result  are  shown  to  help  evaluate  the  grids.  Part  I  of  this  report  is  essentially 
a  self  contained  technical  report  that  is  being  readied  for  submission  to  a  Journal. 

Part  II  of  this  report  is  both  brief  and  sketchy  in  its  presentation.  In  this  section 
we  describe  some  of  our  success  in  treating  bodies  with  sharp  edges  and  bodies 
that  are  exceptionally  concave.  The  hyperbolic  grid  generation  tends  to  propogate 
discontinuities  and  some  bodies  are  not  compatible  with  the  imposed  orthogonality 
and  volume  constraints  that  are  user  imposed.  When  this  happens,  the  hyperbolic 
grid  generator  breaks  down.  What  must  be  done  in  this  case  is  to  relax  these 
constraints.  We  have  had  some  success  doing  this,  but  this  process  remains  still  an 
art  form. 

The  last  part  of  this  report  describes  a  flow  field  algorithm  development.  During 
the  course  of  this  research  we  had  some  considerable  interaction  with  AFWAL, 
and  at  one  point  became  ’side-tracked'  into  a  successful  approach  of  improving  the 
efficiency  of  our  general  implicit  Euler  and  Navier-Stokes  code.  Part  III  of  this 
paper  contains  a  preliminary  paper  describing  this  development. 
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PART  I.  GENERATION  OF 

THREE  DIMENSIONAL  BODY  FITTED  COORDINATES 
USING  HYPERBOLIC  PARTIAL  DIFFERENTIAL  EQUATIONS 


INTRODUCTION 

Body  conforming  curvilinear  grids  are  often  used  in  finite  difference  flow  field 
simulations.  One  reason  for  this  is  that  the  application  of  boundary  conditions 
can  be  simplified  in  finite  difference  calculations  because  grid  lines  coincide  with 
boundary  lines.  This  is  especially  important  in  high  Reynolds  number  viscous  flow 
simulation  in  which  high  flow  gradients  near  the  body  surface  must  be  resolved. 

The  task  of  generating  a  satisfactory  body  conforming  coordinate  system  is  not 
easy.  The  grids  must  not  be  too  distorted,  they  should  have  smooth  variation, 
and  they  should  be  clustered  to  flow  field  action  regions  -  typically  near  boundary 
surfaces.  Moreover,  the  grids  should  be  generated  in  an  automatic  manner  that 
requires  a  minimum  of  user  input. 

One  approach  for  generating  body  conforming  grids  with  minimum  user  input 
has  been  to  solve  a  set  of  partial  differential  equations.  In  this  technique  level  lines 
of  £(*,y,z),  r/(x,  y,  z),  and  f(z,y,  2)  that  have  monotone  variation  are  sought  as  a 
solution  of  a  set  of  partial  differential  equations.  Generally  values  of  £,  9  and  ? 
are  user  specified  on  the  boundary  surface  and  constraints  expressed  as  differential 
equations  are  used  to  develop  the  grid  away  from  the  boundaries.  The  most  pop¬ 
ular  such  approach  requires  the  solution  of  a  set  of  elliptic  equations  that  satisfy 
the  maximum  principle  [l-5j,  however,  hyperbolic  [6, 7)  and  parabolic  [8]  governing 
equations  have  been  used  as  well,  at  least  in  two  dimensional  applications. 

In  this  report  one  way  of  extending  the  hyperbolic  grid  generation  method  of 
Steger  and  Chaussee  [6]  to  three  dimensions  is  developed.  In  two  dimensions  the 
two  differential  constraints 

£*•1*  +  (yQy  =  0  ((fa)) 

h  -  tvv*  =  (avt  1  ((it)) 

or  in  (,  tj  computational  space 

*t*n  +  =  0  ((2a)) 

*{!&»  -  =  AV  ((2b)) 

have  been  solved  by  marching  in  17  from  an  initial  data  plane  9(2, y)  =  constant. 
The  first  equation  is  a  constraint  of  orthogonality.  The  second  equation  controls 
the  mesh  spacing  with  the  user  specifying  the  mesh  control  volume  AV  (actually 
area  in  two  dimensions).  A  linearized  version  of  equations  (2)  is  readily  shown  to  be 


hyperbolic  and  suitable  for  marching  in  9.  Equations  (2)  are  solved  in  computational 
space  to  give  the  x,y  location  of  the  £  =  constant  and  9  =  constant  grid  lines. 

The  two  partial  differential  equations,  expressed  as  either  Equations  (1)  or  Equa¬ 
tions  (2),  have  been  referred  to  as  a  mesh  cell  volume  procedure  for  grid  generation. 
In  the  next  section  a  three  dimensional  extension  of  this  procedure  is  developed. 

THREE  DIMENSIONAL  GRID  GENERATION  EQUATIONS 

A  body  fitted  exterior  grid  about  an  arbitrary  closed  boundary  surface  is  desired. 
Only  a  simple  topology  such  as  that  illustrated  in  Figure  (1)  will  be  considered 
here.  The  body  surface  is  chosen  to  coincide  with  f (z,y,z)  =  0  and  the  surface 
grid  line  distributions  of  £  =  constant  and  9  =  constant  are  user  specified.  The 
outer  boundary  i(x,y,  z )  =  fm42  is  not  specified,  it  is  only  required  to  be  sufficiently 
far  removed  from  the  inner  boundary.  Using  (  as  the  marching  direction,  partial 
differential  equations  are  sought  which  produce  planes  of  constant  £,  9  and  f  to 
form  a  nonsingular  mesh  system. 

An  extension  of  the  mesh  cell  volume  procedure  to  three  dimensions  is  proposed. 
In  three  dimensions,  however,  there  are  three  orthogonality  relations  and  one  cell 
volume  constraint.  At  any  point  four  equations  are  available  to  predict  the  three 
unknowns  x,y  and  z  so  one  equation  must  be  discarded.  Because  (  is  the  marching 
direction  it  is  natural  to  use  only  the  two  orthogonality  relations  that  involve  f,  this 
leads  to  the  governing  equations 

*€*<  +  VtV<  +  =  0  ((3a)) 

+  +  =°  ((3b)) 

x(Vnz<  +  xtVi *n  +  xvV< ~  xtV<zn  ~  xnVt*<  “  x<Vn*t  -  ((3c)) 

or  with  f*  defined  as  (x,  y,  z)* 

'•r‘=°-  IM-r*-AV 

The  first  two  equations  represent  orthogonality  relations  between  £  and  f  and  be¬ 
tween  9  and  f,  while  the  last  equation  is  the  volume  or  finite  Jacobian  constraint. 

Equations  (3)  comprise  a  system  of  nonlinear  partial  differential  equations  in 
which  x,  y  and  z  are  specified  as  initial  data  at  f  =  0.  As  developed  below,  lin¬ 
earization  and  analysis  of  Equations  (3)  about  a  nearby  known  state  reveals  that 
the  system  is  hyperbolic  with  f  as  the  marching  direction. 

Let  x°,y°,z°  represent  a  nearby  known  state  so  that 

*  =  *°  +  (*  -  *°)  =  *°  +  x 


where  x,  y  and  z  are  small.  Substitution  of  these  expressions  into  Equations  (3)  and 
elimination  of  products  of  tilde  terms  results  in  the  locally  linearized  system 


Ao(f-  fo)$  +  Bo(r  -  fo),  4-  Co{f-  r0){  =  f 


with 


A  = 


B  = 


C  = 


xt  ut  z{ 

0  0  0 

itfaz<  ~  V(2v)  {xizn  ~  xvz<)  ixtiV<  ~  x<Vit) 
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0 

v< 


^(v<2€  -V(z<)  ix(z<  -  x<zi)  [x<Vi  -  HV<) 


x( 


V( 

Vn 


[V(2v  - V**t)  ( xnH  -  H zn)  ix(Vv  ~  X*V() 


or .  er 

df  0{ 

Mil.  /=  I  1  = 


0 

0 


AV  -  A  V0  J  V  AV  -  AVb 


•«  . 

Let  R  =  r  -  r0,  then  (5)  is  rewritten  as 

AqR(  +  B0Rn  +  CqR{  =  f 


((5)) 

((6a)) 

((6b)) 

((6c)) 

((6d)) 

m 


Now  C0  1  exists  unless  (AV)"1  — ►  oo,  which  we  will  not  impose,  so  (7)  can  be 
rewritten  as 

C^AoRf  +  C^BqR„  +  R{  =  Cf'f  ((8)) 

Although  the  verification  is  nontrivial,  Cq1Ao  and  Cq'Bq  are  found  to  be  symmet¬ 
ric  matrices  (this  was  carried  out  by  Dennis  Jesperson  of  the  NASA  Ames  Research 
Center,  who  used  MACSYMA).  The  linearized  system  Equation  (8)  is  therefore 
hyperbolic  and  can  be  marched  with  f  serving  as  the  “time-like”  direction. 

It  can  be  pointed  out  that  an  analysis  was  attempted  for  the  three  orthogonality 
relations  alone.  These  equations,  however,  are  found  to  be  improperly  posed  for 
marching  with  initial  data  in  j.  Indeed,  as  best  as  we  can  discern,  the  three  rela¬ 
tions  do  not  lend  themselves  to  unique  solutions  regardless  of  the  type  of  boundary 
conditions  specified. 


**u  v  . 


SOLUTION  PROCEDURE 

The  nonlinear  system  of  grid  generation  equations  given  by  Equations  (3)  are 
solved  with  a  noniterative  implicit  finite  difference  scheme.  An  unconditionally 
stable  implicit  scheme  is  chosen  so  the  marching  step  size  in  f  can  be  arbitrarily 
selected  based  only  on  considerations  of  accurately  generating  the  grid.  Iterative 
solution  of  the  nonlinear  grid  generation  equations  is  avoided  by  expanding  the 
equations  about  the  previous  marching  step.  As  a  consequence  Equation  (7)  is 
solved  with  the  nearby  known  state  0  taken  from  the  previous  f  step, 
a)  Numerical  Method 

Let  A{  =  Aq  =  Af  =  1  such  that  (  =  and  f  =  /  -  1.  Central 

spatial  differencing  of  Equations  (5)  in  (  and  17  with  first  order  backward  implicit 
differencing  in  (  leads  to 

At6((?t+i  —  ?i)  +  £(£,(?<+!  -  fj)  +  CjVfri+i  =  th+i  ((9)) 


where 


{fi+i  = 


" ) 
AVi+1 ) 


S(fj  m 


l  “  *J-1  c  -  *Wl  ~  flk-1 

"j1  4*r*- — 2 

Vffi+i  =  n+i  -  n 


Note  that  CiSffi  was  subtracted  from  /i+i  to  produce  $+1  in  the  above.  Throughout 
only  those  indices  that  change  are  indicated,  thus  ri+j  =>  ry,*,i+i,ry+i  =>  fj+i.k.t 
etc. 

Multiplying  through  by  Cfl  gives 

OflAt8t{n |+i  -  fi)  +  -  r|)  +  I(Fi+i  -  fj)  =  ((10)) 

where  J  is  the  identifying  matrix.  To  reduce  the  inversion  cost  the  difference  equa¬ 
tions  are  approximately  factored  as 


(/  +  C'-'AtScHI  +  CrlBl6„)(r(+l  -r{)  = 


((H)) 


so  that  r}+i  is  obtained  by  solving  sequences  of  one-dimensional-like  block  tridiag¬ 
onal  systems 

{I  +  C,"lA|fi()0*  1+1  =  $1+1  ((12a)) 


(/  +  C^lBfSn)VtFt+i  =  f  <+1 


((12b)) 


m 

-j 


•  A  •  *•*•*•'  •  •  •  •  • •  ■  •  • 


4 


n+i  =  n  +  v{rJ+i 


U12c)) 


In  practise  numerical  dissipation  terms  are  added  in  the  £  and  t)  directions.  Typi¬ 
cally  we  have  used  a  combination  of  fourth  and  second  differences  which  axe  explic¬ 
itly  and  implicitly  included  into  the  basic  algorithm  as 

[/  +  Cf'A,6(  -  £i(AV){]|/  +  Cr'B,i,  -  f,(AV),](n+1  -  F,)  = 

-|£.(AV)|+£.(AV);]F, 

where  tt  is  order  jg  and  e,-  is  2.5  or  more  larger  than  ee.  As  an  alternative  to  fourth 
differences,  second  order  terms  such  as 

e|AVn|AVf$ 

have  been  used  both  implicitly  and  explicitly. 

The  coefficient  matrices  Ai,Bi  and  C\  contain  £  and  i j  derivatives  which  are 
formed  using  central  differences.  These  matrices  also  contain  derivatives  for  x(,y{ 
and  zt  which  are  obtained  from  Equations  (3)  in  terms  of  £  and  t)  derivatives.  That 
is,  Equations  (3)  are  linear  in  the  unknowns  x(,y(  and  z{.  They  are  easily  solved 
for  as 

(*<\  AV  f  ***'  ~  ^ 

*  =  ** " *«*  ((13)) 

\  *C  V  ’  *r,Vt  } 


Det(C)  =  (i/(z,  -  y„z()i  +  (x„z(  -  x(z„)3  +  (x(yn  -  xny()7 

Note  that  AV/^(z*  +  y*  +  z*)  =  Det(C)  so  that  Det(C)  will  be  zero  if  and  only 

If  the  user  specified  AV  =  0.  Hence,  C"1  will  exist, 
b)  Cell  Volume  Specification 

The  user  has  control  of  the  grid  by  means  of  the  initial  surface  point  distribution 
and  by  specification  of  the  cell  volumes,  AVj,*,i .  Through  the  cell  volumes  the 
extent  and  clustering  of  the  grid  can  be  essentially  controlled.  Because  the  cell 
volume  at  each  point  must  be  specified,  it  is  clear  that  the  user  must  devise  some 
kind  of  method  for  determining  volumes.  There  are  many  possibilities,  here  one 
such  approach  is  illustrated. 

Suppose  we  had  a  sphere  to  grid.  A  reasonable  grid  might  have  uniform  angle 
spacing  and  have  a  radial  grid  distribution  that  is  exponential.  For  this  special 
geometry  and  grid  we  can  analytically  determine  the  control  volumes  by  a  simple 
formula.  Now  take  the  problem  at  hand,  perhaps  an  aircraft  fuselage,  which  we 
want  to  mesh  as  a  warped  spherical-like  grid.  We  can  find  a  sphere  that  has  the 
same  surface  area  as  our  fuselage  and  use  the  grid  cell  volumes  of  the  sphere  to 


specify  the  cell  volumes  of  the  fuselage  grid.  However,  the  fuselage  will  not  have 
the  same  kind  of  surface  area  distribution  as  a  sphere  with  equal  angle  distribution. 
So  here  we  need  an  adjustment,  something  like 


—  (A  Vj'iC'i)  sphere  £(l  £) 


ggfcgrg 

(AAy 

,k)jutelage 


((14)) 


where  6  1  for  small  l  and  6  —*  0  for  large  l.  That  is,  the  volumes  would  be 

adjusted  initially  to  the  local  boundary  surface  increments.  But  as  we  march  out 
the  uniform  spherical  volumes  would  gradually  be  specified.  Such  an  approach  has 
been  used,  and,  as  a  result,  the  far  field  portion  of  the  grid  tends  to  be  uniformly 
spherical.  The  results  shown  later  will  illustrate  this  behavior. 

C)  Axis  Treatment 

A  coordinate  singularity  will  be  encountered  whenever  a  regular  grid  is  mapped 
over  a  closed  body  such  that  f  =  0  corresponds  to  the  body  surface.  Here  we 
will  generate  warped  spherical  grids,  so  equation  (5)  becomes  singular  at  the  axis. 
Using  (  as  the  marching  direction  with  £  and  r\  as  sketched  in  Fig.2  (  i.e.  £  from 
pole  to  pole  and  tj  equatorial  )  then  the  axis  at  £  =  0  and  £  =  £mo*  represent 
singularities.  In  particular,  t)  derivatives  approach  zero  at  the  pole  and  equations 
(2b)  and  (2c)  are  lost.  If  there  are  KMAX  points  in  the  ^-direction,  however,  the 
difference  equation  corresponding  to  Eq.(2a)  can  be  imposed  KMAX  times  to  solve 
for  three  axis  unknowns,  namely  x,y,  and  z.  Thus  the  axis  points  are  overdetermined 
and  can  be  solved  for  via  a  generalized  inverse  scheme.  Such  an  approach  has  been 
coded.  However,  it  is  difficult  to  implement  implicitly,  and  it  has  not  proved  to  be 
reliable  for  severely  distorded  axis  cases,  for  example,  when  the  axis  extends  from 
a  wing  tip. 

Rather  than  try  to  compute  the  axis  as  the  solution  evolves,  we  have  instead 
prespecified  the  axis  position.  For  now  a  fixed  straight  axis  has  been  used  that  is 
aligned  with  a  coordinate.  Typically  a  y-coordinate  line  has  been  chosen  so  that  x 
and  z  are  zero  along  the  axis.  Values  of  y  along  the  axis  are  obtained  as  part  of  the 
solution  process  by  imposing  the  orthogonality  conditon,  which  in  this  degenerate 
case  is  simply  yn  —  0.  This  condition  is  imposed  using  second  order  accurate  three 
point  differencing,  and  it  is  implicitly  implemented  into  the  q-block  tridiagonal 
solution  process  using  a  slightly  modified  block  solver. 

RESULTS 

To  demonstrate  the  hyperbolic  grid  solver  a  series  of  grids  were  generated  about 
ellipsoidal  and  'super  ellipsoidal’  body  shapes.  The  body  generating  function  is 
given  by 

plan  form 

(*/*m«*)*  +  [y/Vma*)n  —  I 


thickness 

(t/ j Umax)"1  +  ( z( ^mat)m  —  1 

profile 

(*/*ni4j)  "+■  =  1 

where  n,  m  and  l  are  2  for  ellipsoids  and  4  or  6  or  8  etc.  for  super  ellipsoids.  By 
chosing  very  high  aspect  ratios  the  ellipsoid  can  mimick  a  practical  wing  planform, 
or,  for  more  moderate  aspect  ratios,  a  fuselage.  In  either  case  an  axis  singularity 
exists. 

The  right  half  of  an  elliptic  wing  planform  with  the  user  specified  body  surface 
grid  point  distribution  is  shown  in  Fig.2.  This  'wing'  has  a  16:1  planform  aspect 
ratio  with  a  20  per  cent  thick  ellipse  serving  as  the  airfoil.  Thus  the  body  is  an 
ellipsoid  with  ratios  16  :  1  :  1/5.  A  uniform  grid  spacing  of  1/2  per  cent  thick 
maximum  chord  was  specified  as  the  first  grid  spacing  off  the  body.  Figures.4  to  6 
show  various  segments  of  the  generated  grid  for  £  and  q  =  constant  planes.  Figure 
4  shows  an  q  =constant  plane  (  here  in  the  plane  of  the  wing  trailing  edge)  at  the 
wing  tip.  As  seen,  the  grid  is  very  smooth,  uniform  off  the  body,  and  the  axis  at 
*  =  0  is  well  behaved.  A  lower  half  of  the  grid  at  the  wing  midspan  (  y  =  0)  is  shown 
in  Fig.  5.  This  is  a  £  =  constant  plane.  Again  the  smoothness  and  grid  clustering 
control  is  illustrated.  Also  in  the  far  field  the  grid  is  tending  to  be  spherical  because 
spherical  mesh  incremental  volumes  have  been  specified.  Finally,  views  above  the 
wing  in  the  x  =  0  plane  are  shown  in  Figs.  6a  and  6b.  These  views  are  focused  at 
the  wing  tip  and  again  illustrate  satisfactory  axis  treatment. 

A  similar  grid  was  computed  as  illustrated  by  Figs.  7  to  10.  The  wing  in  this 
case  used  a  more  realistic  airfoil  shape  that  has  a  thickness  ratio  of  15  per  cent. 
The  grid  views  Figs.  7  and  9  give  an  indication  of  the  wing  and  airfoil  section.  Note 
that  the  airfoil  has  a  rounded  trailing  edge. 

It  must  be  remarked  that  if  the  airfoil  trailing  edge  radius  is  continuously  reduced 
that  the  method  breaks  down.  "Pretty  bows"  get  tied  within  the  grid.  Provided 
the  trailing  edge  radius  is  not  zero,  break  down  can  usually  be  avoided  by  clustering 
grid  points  in  this  region  and  adjusting  the  specified  volumes.  In  lieu  of  good  surface 
clustering  functions  one  can  sometimes  obtain  an  adequately  resolved  trailing  edge 
by  generating  an  overall  much  finer  grid  than  what  is  desired  for  the  flow  solver. 
In  this  case  we  simply  discard,  say,  every  other  grid  point  to  obtain  the  final  grid. 
Because  the  hyperbolic  grid  solver  is  quite  efficient,  we  can  readily  operate  in  this 
manner.  Ultimately  the  scheme  does  break  down  when  subject  to  a  truely  sharp 
trailing  edge. 

The  above  grids  have  grid  spacing  that  is  adequate  for  inviscid  flow  simulations. 
The  hyperbolic  grid  generation  procedure  can  also  directly  generate  grids  suitable 
for  high  Reynolds  number  viscous  flow  simulation  by  simply  having  the  user  specify 
a  much  finer  mesh  spacing  going  out  away  from  the  surface.  In  Figs.  11  and  14  we 


show  views  of  a  grid  generated  about  a  bluff  body  fuslage.  The  body  in  this  case 
is  defined  as  a  2  :  1  ellipsoid  in  the  x-y  and  y-z  planes  (  see  Fig.  11).  In  the  x-z 
plane  a  1 :  1  super  ellipsoidal  cross  section  in  specified  by  seting  m  =  4  (  see  Fig. 
14).  The  first  grid  spacing  off  the  body  in  this  case  was  specified  as  0.00002  of  the 
body  cross  sectional  diameter.  Magnified  grid  lines  coming  into  the  axis  are  shown 
in  Fig.  13  to  indicate  this  fine  grid  spacing  just  off  the  body. 

Although  a  less  interesting  configuration,  a  similar  viscous  grid  was  generated 
about  a  6  :  1 : 1  ellipsoid.  To  demonstrate  that  the  grids  being  generated  with  this 
approach  are  indeed  useful,  Figs.  15  to  16  show  calculated  particle  paths  and  surface 
’oil  flow’  simulations  on  this  6:1:1  ellipsoid.  These  Naiver  Stokes  calculations 
were  carried  out  by  T.H.  Pulliam  using  ARC3D.  The  free  stream  conditions  were 
chosen  as  Mach  number  of  0.74,  angle  of  attack  of  25  deg.,  and  0  yaw  angle  with 
a  plane  of  symmetry  imposed.  The  Reynolds  number  was  44  x  106  based  on  the 
diameter  and  laminar  flow  was  assumed. 


CONCLUSION 

A  procedure  has  been  developed  for  generating  body  fitted  coordinates  in  three 
dimensions  using  hyperbolic  partial  differentia]  equations.  In  this  report  the  scheme 
has  been  used  to  generate  warped  spherical-like  grids  about  simple  ellipsoidal  wing 
and  fuslage  configurations.  The  hyperbolic  grid  generator  can  fail  whenever  the 
body  surface  is  discontinuous  or  the  user  specified  surface  grid  distribution  is  too 
irregular.  For  simple  continous  body  shapes,  however,  the  hyperbolic  grid  genera¬ 
tion  scheme  can  be  very  fast  and  reliable.  It  requires  a  minimum  of  user  interaction, 
and  it  can  be  used  to  generate  grids  suitable  for  either  inviscid  or  viscous  flow  sim¬ 
ulation. 
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Figure  1.  Well-Ordered  Warped  Spherical  Grid  Mapping 
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PART  II.  EXTENSIONS  TO  HYPERBOLIC  GRID  GENERATION 

In  this  very  brief  section  we  shall  simply  make  mention  of  some  problem  areas 
and  some  'fixes’  in  dealing  with  configurations  with  sharp  edges  or  exceptionally 
concave  cross  sectionals.  In  these  situations  the  hyperbolic  grid  generation  scheme 
generally  breaks  down.  Figure  1  indicates  what  breakdown  can  look  like.  Only  two 
dimensional  results  are  considered  here  with  (  the  coordinate  along  the  body  and 
q  the  marching  coordinate. 

For  exceptionally  concave  bodies  T.  J.  Barth  found  that  he  could  often  make  the 
hyperbolic  grid  generator  still  perform  by  multiplying  the  ^-difference  by  a  factor 
greater  than  one.  This  ad  hoc  fix  seems  to  allow  the  generated  grid  lines  more 
flexibility  to  a4just  in  £,  yet  does  not  seem  to  too  adversely  effect  orthogonality. 
The  success  of  this  approach  is  indicated  by  Figures  2  to  4  which  show  cross  sectional 
grid  views  about  the  X-24c  which  where  generated  for  AFWAL.  Barth  and  Risk  have 
also  generated  grids  about  cross  sections  of  the  Shuttle,  a  grid  is  shown  in  F>"  5. 

The  hyperbolic  grid  generation  scheme  will  also  break  down  when  a  grid  must  be 
wraped  around  a  sharp  edges  such  as  the  trailing  edge  of  an  airfoil  meshed  with  an 
"O-grid”  (Le.  warped  cylindrical  grid).  However  by  treating  the  grid  line  leaving 
the  trailing  edge  as  a  special  surface  at  which  orthogonality  is  not  imposed,  airfoil  O- 
grids  have  been  generated.  Figures  6a  and  6b  show  a  grid  about  a  cambered  airfoil 
is  which  the  (=constant  line  emanating  from  the  trailing  edge  is  specially  treated. 
Near  the  body  this  line  is  determined  so  that  it  has  a  slope  that  bisects  the  trailing 
edge  angle,  while  the  arc  length  increments  between  points  is  user  specified.  Away 
from  the  body  surface  this  special  logic  is  blended  into  the  usual  grid  generation 
formula.  As  the  grid  views  indicate,  such  trailing  edge  logic  can  work  well,  but 
proper  specification  of  the  arc  lengths  along  the  special  plane  is  still  very  much  an 
art.  We  are  also  working  on  special  logic  that  uses  slope  and  volume  information, 
and  this  is  evolving  into  a  more  automatic  process.  A  result  obtained  by  Risk  and 
Barth  is  shown  in  Fig.7.  Some  irregularity  can  still  be  obtained,  but  this  process  is 
becoming  reliable  and  automatic. 
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Figure  5.  Grid  About  Shuttle  Cross  Section 


Figure  7.  Grid  About  Cambered  Airfoil 
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PART  III.  An  Efficient  Approximate  Factorization  Implicit  Scheme 

for  the  Equations  of  Gasdynamics 

Timothy  J.  Barth*  and  Joseph  L.  Steger* 

NASA  Ames  Research  Center,  Moffett  Field ,  CA  94035 


L  Introduction 

Numerical  algorithms  for  solving  the  compressible  Euler  and  Navier-Stokes 
equations  have  received  considerable  attention  over  the  last  decade.  Both  im¬ 
plicit  and  explicit  time  advance  and  iterative  techniques  have  been  ultilized. 
Explicit  methods  typically  have  lower  arithmetic  operation  counts  but  more 
restrictive  stability  bounds  than  implicit  methods.  Use  of  implicit  methods  is 
generally  advantageous  for  flow  situations  where  an  explicit  time  step  limita¬ 
tion  would  be  much  less  than  the  time  step  neccessary  for  the  desired  ac¬ 
curacy.  This  situation  frequently  arises  in  the  solution  of  unsteady  flows 
that  are  characterized  by  low  frequency  motion,  boundary  condition  forcing, 
and  high  Reynolds  number  viscous  effects.  Implicit  algorithms  can  often  be 
readily  adapted  to  be  efficient  relaxation  schemes  for  steady  state  applications 
as  well. 

The  purpose  of  this  paper  is  to  present  a  technique  that  substantially 
reduces  the  arithmetic  operation  count,  but  otherwise  retains  the  stability 
characteristics  of  a  Beam- Warming  approximate  factorization  implicit  al¬ 
gorithm  for  the  Euler  equations.  This  new  method  extends  a  matrix  splitting 
of  Steger[l],  previously  restricted  to  Cartesian  coordinates,  such  that  equa¬ 
tion  decoupling  is  achieved  for  the  conservative  form  of  the  Euler  equations. 
The  equation  decoupling  (  or  matrix  reducibility  )  results  in  a  significant  im¬ 
provement  in  the  efficiency  of  the  algorithm.  By  ultilizing  a  local  transforma¬ 
tion,  the  equation  decoupling  techniques  are  extended  to  the  Euler  equations 
in  generalized  coordinates.  Computation  of  transonic  flow  about  a  NACA 
0012  airfoil  is  used  to  verity  that  no  numerical  stability  is  lost  with  the  new 
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reducible  implicit  algorithm.  Extension  of  the  technique  to  the  Navier-Stokes 
equations  is  also  indicated  in  the  Appendix  A. 


IL  Background 

In  this  section  we  review  a  sound  speed  and  velocity  splitting  concept 
for  an  approximate  factorization  implicit  algorithm  in  Cartesian  coordinates 
and  describe  how  it  can  be  used  to  improve  computational  efficiency.  The 
extension  for  general  coordinate  systems  is  described  in  the  following  section. 

a)  Beam  and  Warming  Algorithm 

The  conservative  form  of  the  Euler  equations  for  a  perfect  gas  can  be 
written  in  Cartesian  coordinates  as 


dtQ  "1"  "I"  dy F  —  0 
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Here  p  is  the  fluid  density,  u  and  v  are  the  Cartesian  velocity  components,  e 
is  the  total  energy  per  volume,  and  p  is  the  pressure. 

A  general  purpose  implicit  finite  difference  scheme  (c.f.  Refs.  2  and  3)  for 
equation  (2.1)  that  can  be  readily  adapted  to  either  steady  or  unsteady  flow 
applications  is  given  by 


I  +  h6xAn  -h  I  +  h6yBn  +  AQn  -- 


-At 


(2.4) 


where 


Jj(4)  =  e,At[(VAg  +  (VA^] 


(2.5a) 


and 

D&  =  -etAf(VA)„  D^  =  -f»A«(VA)y  (2.5b) 

Here  5,  and  5y  are  three-point,  second  order  accurate  central  space  difference 
operators,  while  A  and  V  denote  forward  and  backward  space  differences. 
Either  first  or  second  order  accurate  time  differencing  may  be  used  with  h=At 
or  At/2,  alternately  the  time  step  may  be  used  as  a  relaxation  parameter 
and  may  vary  in  space  as  well.  With  the  use  of  central  space  differencing 
numerical  dissipation  is  needed  and  implemented  with  the  second  and  fourth 
order  differences  given  above.  The  parameters  ee  and  e*  control  the  amount  of 
numerical  dissipation.  These  parameters  may  be  constants  or  may  vary  with 
the  solution  gradients  (  in  which  case  they  are  moved  inside  the  operators  so 
as  to  maintain  conservative  form).  Steady  state  convergence  can  be  greatly 
enhanced  by  more  implicit  treatment  of  the  dissipation  terms  and  by  use  of 
space  varying  scaling  parameters. 

In  equation  (2.4),  local  time  linearizations  have  been  ultilized  to  avoid 
solving  nonlinear  equations  between  time  or  iterative  levels  n  to  n-f-1  by 
expanding  E  and  F  in  terms  of  A Q  as 

£n+1  =  En  -+-  An(Qn+1  -  Qn ),  A  =  dE/dQ  (2.6a) 

Fn+i  __  Fn  +  sn(Qn+l  -  Qn),  B  =  dFfdQ  (2.6b) 

The  left  hand  side  operators  of  equation  (2.4)  form  a  linear  system  of 
equations  that  must  be  solved  at  each  time  or  iteration  level.  Although 
the  matrix  band  structure  of  this  system  has  been  altered  by  approximately 
factoring  into  one-dimensional-like  operators,  the  inversion  process  has  been 
simplified.  The  approximate  factorization  (AF)  reduces  the  inversion  work 
but  sometimes  at  the  cost  of  limiting  the  time  step  At  because  of  either 
time  accuracy  considerations  or  because  of  reduced  inefficiency  in  relaxation 
applications. 

In  practice  the  difference  equations  (2.4)  are  solved  in  the  ADI  algorithm 
form 

L*AQ*  =  RHS(2A)  (2.7a) 

IyAQ"  =  A  Q*  (2.7b) 

Variants  of  this  algorithm  include  higher  order  difference  approximations 
(including  pseudo-spectral  right  hand  side-operators  (4)  ),  space  varing  At 
for  steady  state  applications,  addition  of  viscous  terms,  etc  [c.f.  5-9  ]. 


b)  Diagonalization 

Efforts  have  been  made  to  reduce  the  inversion  work  over  and  above 
what  is  obtained  with  approximate  factorization.  In  particular,  Pulliam  and 
Chaussee7  have  used  eigenvector-similarity  transforms  to  reduce  the  block 
tridiagonal  matrices  to  scalar  tridiagonals(  see  also  Chaussee  and  Pulliam8, 
Steger,  Pulliam  and  Chima8,  Coakley  10  ).  Specifically,  the  eigenvectors  of  the 
A  and  B  Jacobian  matrices,  denoted  as  X  and  Y ,  are  further  approximately 
factored  into  the  left  hand  side  of  equation  (2.4)  as 
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(2.8) 


where 

Aa  =X~lAX  and  A  b  =  Y~1BY 

As  Figure  1  illustrates,  the  number  of  operations  in  solving  a  block  tridiagonal 
(  that  uses  only  L-U  decomposition  )  increases  with  the  cube  of  the  block 
size,  a  formula  given  by  Merriam  11  indicates  the  work  increases  as  ^m3  -f- 
| m2—  where  m  is  the  dimension  of  the  block.  Specifically  our  own  coded 
4x4  block  solver  requires  370  operations  per  entry  (i.e.  grid  point)  while  the 
lXl  tridiagonal  needed  in  (2.8)  requires  only  9  operations  per  entry  (or  36 
operations  for  4  scalar  tridiagonals).  Although  the  eigenvector  matrices  must 
be  formed  and  multiplied  through  (each  4x4  multiply  requires  28  operations 
per  entry),  the  arithmetic  operation  count  of  equations  (2.8)  is  considerably 
reduced  from  that  of  equations(2.4).  For  block  tridiagonal  matrices  subject 
to  periodicity  conditions  and  for  block  pentadiagonal  matrices  the  saving  is 
even  more  significant. 

e)  Sound  Speed  and  Velocity  Splitting 

An  alternate  method  for  reducing  the  inversion  (i.e.  solution  )  work  was 
proposed  and  demonstrated  in  reference  [l].This  method  was  also  motivated 
by  similarity  transforms,  but  it  does  not  explicitly  use  them  in  the  algorithm. 
Instead,  it  involves  splitting  the  Jacobian  matrices  A  and  B  into  velocity  and 
sound  speed  (  or  pressure  parts  )  as 

A  =  An  “l-  At 
B  =  Bv  -f*  Bc 


(2.9a) 

(2.9b) 


such  that  the  eigenvalues  a(A)  and  o(B)  are 


with 


<j{A)  =  a{Au)  +  a{Ac) 
ff[B)  =  ff{Bu)  +  o{Bc) 

<j(Au)  =  ( u ,  u,  u,  u),  a(Ac)  =  (0,  c,0,  —  c) 
o(Bv)  =  (v,  v,  v,  v),  a(Bc)  =  (0, 0,  c,  —  c) 


(2.10a) 

(2.10b) 

(2.11a) 

(2.11b) 


The  matrices  Au,  Ac,  Bv  and  Bc  were  found  from  a  natural  reduced  form  of 
the  equations  in  nonconservative  from.  Specifically  Ac  and  Bc  were  split,  as 


Ac  =  (7—1) 
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where 


a4i  =  W(«2  +  «2)/2]  -  7«P/Ip(7  -  l)2] 
=  IP  Ml  ~  I)2!  “  u2 
b\x  =  [v(u2  +  v2)/2]  -  7 ~  l)2] 
643  =  7P/ W7  -  l)2]  “  1,2 


while  A«  and  B«  are 


A^u  —  A.  ^ 
Bv  —  B  —  Bg 


(2.13a) 

(2.13b) 


This  splitting  produces  matrices  A«  and  Bv  that  are  more  complex  than  A 
and  B.  However,  Steger  [1]  noted  that  Q  is  an  eigenvector  of  both  matrices, 
i.e. 


AUQ  =  uQ 
BVQ  =  vQ 

which  prompted  the  ad  hoc  substitution 


(2.14a) 

(2.14b) 


V 

A«=(u/+Ac)Q 

(2.15a) 

—  —4 

m 

BQ  =  ( v I  +  BJQ 

(2.15b) 

The  matrices  ul  Ac  and  vl  -}-  Bc  have  a  reduced  form  that  simplifies 
inversion  compared  to  A  and  B. 

Insertion  of  equation  (2.15)  into  the  equations  for  local  linearization  of  the 
Jacobians  (2.6)  produces 

% 

En+\  =  En  _j_  (w/  _j_  Ac)n(Qn+l  -  Qn)  (2.16a) 

pn+l  =  F n  _J_  (y/  _J_  £c)n(Qn  +  1  _  Q")  (2.16b) 

XJltilizing  these  linearizations  in  the  basic  algorithm  equation(2.4) 

LxLy&Q  =  RHS(2A)  (2.17) 

gives  new  left  hand  side  operators  Lx  and  Ly  that  are  easier  to  solve 


Lx  = 


I  +  M,(u/  +  Ac)"  +  £><?>] 

l,=  i + hsy(„r + Be)" + 


(2.18a) 

(2.18b) 


The  end  result  of  this  splitting  is  that  the  new  operators  Lx  and  Ly  form 
matrices  that  no  longer  require  4x4  block  tridiagonal  inversions.  Leaving 
out  the  dissipation  terms  to  illustrate  the  structure  of  these  operators,  we 
obtain 
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where  ae  and  be  are  the  respective  elements  of  Ac  and  Bc  given  by  equation 

(2.12). 

For  the  Lx  operator,  for  example,  the  first  and  third  rows  decouple  from  the 
system  and  can  be  solved  as  scalar  tridiagonal  matrices  with  their  respective 
right  hand  sides.  Once  these  rows  are  solved,  the  elements  of  the  first  and 
third  columns  can  be  moved  to  the  right  hand  side.  The  second  and  fourth 
equation  remain  coupled  and  are  solved  as  a  2  x  2  block  tridiagonal  matrix. 
The  dissipation  terms  also  form  diagonal  tridiagonals  and  so  do  not  alter  the 
structure.  The  block  decoupling  of  the  Ly  operator  is  even  more  conspicuous 
and  is  inverted  (i.e.  solved  for  )  in  a  similar  manner. 

Substitution  of  the  left  hand  side  operators  given  by  (2.18)  in  place  of 
those  of  (2.7)  results  in  a  substantial  reduction  in  arithmetic  operations.  Our 
typical  2x2  block  tridiagonal  requires  55  operations  per  point,  so  the  overall 
inversion,  including  the  two  scalar  tridiagonals,  requires  73  operations  per 
entry.  Because  the  two  scalar  tridiagonals  have  identical  coefficients  this  work 
can  be  even  further  cut  by  solving  them  together. 

The  matrix  splitting  (2.15)  produces  the  flux  vectors 


E  =  AQ  —  uIQ  -f*  ACQ  —  Ev  -f-  Ec 
F  =»  BQ  =  vIQ  +  BeQ  =  FV  +  Fe 


(2.20a) 

(2.20b) 


where 
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(2.21) 


Note  that  we  do  not  reproduce  the  matrix  splitting  (2.9)  (  with  Ac  and  Be 
defined  from  (2.12) )  by  taking  the  Jacobians  of  (2.21).  Surprisingly,  linear 
stability  analysis  (  see  Appendix  B)  as  well  as  numerical  experiment  have 
shown  that  the  use  of  the  Jacobian  matrices  for  Ac  and  Bc  is  unsatisfac¬ 
tory.  However,  the  same  stability  analysis  and  numerical  experimentation 
have  confirmed  that  no  numerical  stability  degradation  occurs  by  using  the 
substitute  linearization  matrices  (2.15).  But,  the  linearization  (2.16)  is  only 
first  order  accurate  because  of  the  substitution  (2.14). 

The  matrix  splitting  concept  can  be  extended  to  three  dimensions.  In  this 
case  the  difference  equations  can  be  represented  by 


LxLvL,{Qn+l  -  Qn)  =  - &t[6xEn  -j-  SyFn  -f-  8,Gn  +  D^Qn J  (2.21) 


Associated  with  each  operator  are  5  X  5  block  tridiagonal  matrices,  with  each 
matrix  inversion  requiring  about  700  operations  per  grid  point.  Typically  the 
inversion  represents  70  to  80  per  cent  of  the  overall  work  per  point.  Replacing 
each  Jacobian  A  by  ul  -j-  Ac  etc.  reduces  the  inversion  cost  from  700  to 
82  operations  per  entry  of  each  block  tridiagonal.  By  combining  the  scalar 
tridiagonals,  this  can  be  cut  to  74  operations. 


HI.  Sound  Speed  and  Velocity  Splitting  in  Steady  Generalized  Coordinates 

The  two  -  dimensional  Euler  equations  can  be  transformed  from  Cartesian 
coordinates  to  steady  curvilinear  coordinates  using  new  independent  variables 

r  =  t 

t  =  Z{x,y)  (3.i) 

v  =  v{z,y)- 


The  Euler  equations  written  in  steady  generalized  curvilinear  coordinates 
are 

dTQ  +  d^E  +  dnF  =  0  (3.2) 
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where  J  is  the  Jacobian  £xt]y  —  £yijx  and 


U=£xu- \-ZyV,  V  =  tixu-\-v  vv 

are  unsealed  contravariant  velocity  components. 

Application  of  the  AF  implicit  scheme  (2.4)  to  equation  (3.2)  is  given  by 


A  Q  = 


where 

and 


/  +  hS'A*  +  Z^2)][/  +  h5vBn  +  D& 

-  +  6nFn  +  DWQn 

2>(4)  =  feAt/-l{(V^ A^)2  -f-  (V7A,)2J/ 


Df  =  -eAtJ-'VzAtJ,  D&  =  —€iAtJ~~1Vt)AT)J 


(3.3a) 

(3.3b) 

(3.3c) 


The  Jacobian  matrices  of  E  and  F  can  be  given  in  terms  of  the  Jacobian 
matrices  A  and  B  as 


A=e,A+|yB 


(3.4a) 


(3.4b) 


B  =  r\xA  +  rjvB 

Making  the  substitutions  of  (2.15  )  into  AQ  and  BQ  results  in 

AQ  =  ( Ut  +  U-A.  +  f  ,BC)Q  (3.5a) 

BQ  =  (V/  +  l uA<  +  ti,Bc)Q  (3.5b) 

Because  of  the  linear  combination  of  both  Ac  and  Bc  being  present  in 
Eq.(3.5),  only  3x3  reducibility  is  achieved.  When  applied  in  the  implicit 
AF  algorithm,  this  results  in  a  savings  as  indicated  in  Ref.l,  but  not  as  great 
a  saving  as  what  was  achieved  in  Cartesian  coordinates.  Moreover,  in  three 
dimensions  one  can  only  reduce  to  a  4  X  4  block  tridiagonal. 

The  transformed  equations  given  above  have  only  been  transformed  in  their 
independent  variables.  That  is,  Cartesian  velocity  or  momentum  components 
have  been  kept,  and,  as  a  result,  the  so  called  strong  conservation  law  form  is 
maintained.  Suppose,  for  example,  we  had  used  orthogonal  body  coordinates 
s  and  n,  and  that  we  transformed  the  momentum  components  as  well.  We 
would  then  obtain  equations  with  coordinate  source  terms,  but  otherwise 
the  equations  would  look  like  their  Cartesian  counterparts  and  the  Jacobian 
matrices  could  be  reduced  as  before.  We  can  achieve  this  same  effect  with 
equation  (3.2)  by  multiplying  through  by  the  matrix 

1  0  0  0- 

c_  o  u  e,  0 

0  v*  'll  0 
-0  0  0  1- 

This  matrix  transforms  u,v  momentum  components  to  U,V  components. 
With  multiplication  of  (3.2)  by  C  we  obtain 
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with 

H=-c[(cer‘B+(c,r'F| 

These  equations  look  much  like  their  Cartesian  counterparts  except  for 
the  coordinate  source  term  and  the  appearance  of  two  extra  pressure  terms. 
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However,  for  orthogonal  coordinates,  V£  •  Vr]  =  0,  and  the  extra  pressure 
terms  then  drop  out.  Although  we  will  not  give  the  details  here,  one  can 
anticipate  from  the  form  of  the  equations  that  the  flux  Jacobians  in  the 
implicit  algorithm  are  2x  2  reduciable  (  here  the  Jacobian  matrices  are  formed 
with  respect  to  pU  and  pV  components,  not  pu  and  pv  components).  So  if 
the  coordinate  generated  source  term  is  either  treated  explicitly  or  properly 
distributed  in  the  L ^  or  Lv  operators,  the  inversion  efficiency  of  the  previous 
section  can  be  obtained. 

Generally  the  transformed  coordinates  will  not  be  orthogonal  and  multi¬ 
plication  through  by  C  will  not  lead  to  the  2x2  reduced  matrix  operators. 
However,  if  we  use  a  transform  matrix  that  brings  out  a  mixture  of  con- 
travariant  and  covariant  velocity  components,  the  extra  pressure  terms  will 
also  disappear.  Although  it  may  not  be  immediately  obvious,  we  can  again 
obtain  the  reduced  matrix  operator  structure  without  requiring  the  coor¬ 
dinates  to  be  orthogonal. 

Consider  the  transform  matrix  C  given  by 
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where 

h  =  v/vf  *  V£ 
1 2  =  >jvr)  ■  V»7 


(3.8) 


and  the  inner  2  x  2  of  C  transforms  n  and  v  into  the  covariant  velocity 
components 


U  =  -  n,v),  V  =  £T‘(-  fy»  +  f.») 

Since  the  determinant  of  C  is  its  inverse  exists  for  all  mappings  of 

interest  and  is 
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Multipling  equation  (3.2)  by  C  we  obtain 


where 


drQ  +  d^E  +  d„F  =  H 


H  =  -C(Ci  'k  +  c,  'f) 


(3.10) 


and  Q  =  CQ,  E  =  CE,  F  =  CF.  The  flux  vectors  written  out  are 
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(3.11) 


Unlike  equation  (3.7)  extra  pressure  terms  are  not  picked  up  in  the  momentum 
equations. 

The  previous  numerical  algorithm  extended  to  equation  (3.10)  is  given  by 


+  >|/  +  hS„B"  +  Df>  Ag"  = 

-  At[f(£*  +  SnF"  +  Hn  +  D(4>Q' 


(3.12) 


where  we  have  lagged  the  source  term  and  where  A  =  dE/dQ  and  B  = 
OF/dQ.  _  _ 4  _ t 

It  can  be  verified  that  A  —  CAC  and  B  =  CBC  .  Using  these 
relations  as  well  as  (3.5)  we  find 

AQ  =  (UI  +  Ac)Q  BQ  =  (VI  +  Bc)Q  (3.13) 

with 
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and  _  . 

£12  =  \/Vf  •  Vj? 

Making  these  substitutions  into  equation  (3.12)  we  obtain  the  reduced  form 
of  the  algorithm 

/  +  h6£UI  +  Ac)n  +  D^lf/  +  h6„(VI  +  Sc)n  +  AQn  = 

jr  J  •,  (3.15) 

-  At  6tEn  -f  +  Hn  -f-  D&Q 

Finally,  to  avoid  the  source  term  as  well  as  to  take  advantage  of  existing 
codes,  one  can  rewrite  these  equations  in  the  form 

c~l  i + hs€(ui + A*)" + Df> ]  i  +  h6n{vi  +  Bc)n + nW jcAQn  = 

-  Af  + $„f” + 

(3.16) 

The  right  hand  side  of  this  equation  is  the  same  as  equation  (3.3a),  while 
the  left  hand  side  maintains  the  update  for  Q  quantities  rather  than  Q. 
The  advantage  of  this  last  form  is  that  one  can  readily  take  an  existing 
code  for  equation  (3.2)  and  modify  only  the  left  hand  side  operators  so  as  to 
take  advantage  of  the  reduced  inversion  work.  The  steady  state  solution  of 
equation  (3.16)  is  clearly  identical  to  that  of  equation  (3.3a). 


IV.  Results 

The  algorithm  given  by  equation  (3.16)  was  tested  on  a  NACA  0012  airfoil 
under  transonic  flow  conditions  in  order  to  verify  that  no  adverse  stability 
effects  are  incurred  by  using  the  matrix  reductions.  Versions  of  the  basic 
NASA  Ames  Research  Center  AIR2D/ARC2D,  which  solve  equation  (3.3a), 
were  modified  to  the  form  of  equation  (3.16).  As  noted  earlier,  only  that 
portion  of  the  code  which  deals  with  the  left  hand  side  operator  has  to  be 
altered  in  order  to  implement  the  new  algorithm.  The  basic  code  is  described 
in  Refs.  [5,6]. 

A  “C-type"  topology  was  used  for  the  airfoil  calculations.  The  grids  were 
generated  with  either  a  hyperbolic  grid  solver  12>13  that  imposes  orthogonality 
up  to  numerical  errors  of  truncation  and  smoothing,  or  an  algebraic  construc¬ 
tion  that  uses  a  method  similar  to  the  control  function  approach  of  Eiseman14. 
In  all  cases,  the  far  field  boundaries  were  placed  approximately  20  chords 
from  the  body.  Boundary  conditions  on  the  body,  wake,  and  freestream  are 
further  described  in  Refs.[5,6]. 

As  a  first  calculation,  a  relatively  coarse  grid  was  used  which  was  generated 
with  the  hyperbolic  grid  solver.  This  grid  has  157  points  in  the  ^-direction 
and  33  points  in  the  ^-direction  (see  Figures  2a-2b).  The  pressure  distribution 
on  the  NACA  0012  for  M, =  .8  and  a  —  0.  is  shown  in  Figure  2c.  This 
distribution  was  computed  using  both  the  standard  algorithm  (3.3a)  and  the 
reduced  matrix  form  (3.16).  Both  algorithms  were  run  using  Euler  implicit 
differencing  which  is  first  order  accurate  in  time.  As  shown  in  figure  2d,  the 
residual  histories  are  virtually  identical.  However,  the  standard  algorithm 
required  120  CPU  seconds  of  CRAY-XMP  time  per  1000  iterations  while  the 
new  algorithm  required  54  CPU  seconds  per  1000  iterations. 

Numerical  experimentation,  in  which  At  was  adjusted  over  a  wide  range, 
has  confirmed  that  the  new  algorithm  is  as  stable  as  the  standard  algorithm. 
As  an  example  of  the  robustness  of  the  formulation  using  (3.16),  a  much  finer 
grid  was  used  to  compute  the  previous  case.  Figures  3a-3b  show  a  view  of 
the  grid  with  249  points  in  the  ^-direction  and  50  points  in  the  ^-direction. 
The  solution  on  this  grid  (Figures  3c-3d)  was  computed  using  a  spacially 
varying  time  step  scaled  by  -^-j=  where  J  is  the  metric  Jacobian.  For 

this  calculation,  the  drag  coefficient  was  both  converged  to  5  digits  in  600 
iterations  or  68  seconds  of  CRAY-XMP  time.  Figure  3e  shows  the  residual 
history  for  this  calculation. 

Recently,  Beam  and  Bailey  have  reported  in  private  communication  [see 
also  6]  that  significantly  faster  convergence  can  be  obtained  by  incorporat- 


ing  the  fourth  order  smoothing  terms  implicitly  into  the  Beam  and  Warming 
algorithm.  For  the  standard  algorithm,  imposing  fourth  order  dissipation  im¬ 
plicitly  necessitates  4x4  block  pentadiagonal  inversions  which  are  relatively 
expensive.  The  reduced  matrix  algorithm,  however,  requires  only  scalar  and 
2x2  block  pentadiagonal  inversions  which  are  much  less  costly.  In  fact, 
a  special  2x2  block  pentadiagonal  inversion  algorithm,  coded  for  this  al¬ 
gorithm,  requires  only  107  arithmetic  operations  per  entry. 

In  testing  the  sound  speed- velocity  splitting  algorithm  with  implicit  fourth 
order  dissipation,  a  more  advanced  form  of  numerical  filter  was  implemented. 
This  filter  uses  both  second  and  fourth  order  differences  such  that  the  second 
order  difference  is  dominant  near  shock  waves.  The  spectral  radius  of  the 
flux  Jacobian  matrices  is  used  as  a  scaling  to  the  smoothing  coefficients  in  an 
attempt  to  mimick  the  dissipative  nature  of  second  order  flux  split  upwind 
schemes  1S.  Full  details  are  given  in  reference  [6]. 

As  a  test  case,  an  algebraic  grid  was  used  to  compute  the  flow  field  around 
a  NACA  0012  airfoil  at  Mqq  =  .8  and  a  =  1.25  degrees.  The  grid  has 
193  points  in  the  ^-direction  and  33  points  in  the  ^-direction.  The  grid  and 
flow  field  solution  are  shown  in  figures  4a-4d.  The  residual  history  (figure 
4e)  shows  the  rapid  convergence.  For  this  particular  case,  the  CPU  time  on 
the  CRAY-XMP  was  93  seconds  per  1000  iterations.  This  compares  with  65 
seconds  per  1000  iterations  using  the  tridiagonal  version  of  the  new  algorithm 
with  constant  smoothing  coefficients  which  had  much  slower  convergence.  It 
should  be  noted  that  the  smoothing  filter  used  for  the  pentadiagonal  version 
required  16  seconds  more  CPU  time  than  the  same  smoothing  with  a  constant 
coefficient  filter. 

Our  computational  results  verified  that  the  reduced  algorithm  was  every  bit 
as  numerically  stable  as  the  standard  algorithm.  This  experience  is  similar  to 
what  was  observed  previously  by  Steger  on  a  algorithm  version  that  reduced 
to  a  3  X  3  block  and  a  scalar.  Here,  with  reduction  to  the  2  x  2  block,  just 
over  a  factor  of  two  was  saved  in  overall  computer  time  by  using  the  reduced 
algorithm  (3.16)  in  place  of  (3.3a).  The  precise  savings  that  can  be  obtained 
is,  of  course,  machine  and  coding  dependent  and  a  larger  improvement  should 
be  obtained  in  three  dimensions. 

V.  Conclusions 

By  substituting  similar  reducible  matrices  for  the  Jacobian  matrices  of 
the  flux  vectors,  a  substantial  reduction  has  been  achieved  in  the  arith¬ 
metic  operations  needed  in  solving  approximate  factorization  implicit  al- 


t 


gori thins  such  as  Beam-Warming.  Numerical  calculations  indicate  that  the 
new  similarity  split  matrices  can  apparently  be  used  without  any  loss  of 
numerical  stability,  although  time  accuracy  can  be  degraded  unless  additional 
steps  are  taken.  The  new  method  can  be  readily  retrofit  within  existing  codes, 
and  can  likely  be  extended  to  viscous  flow  calculations  as  well. 
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Appendix  A.  Extension  of  Sound  Speed  and  Velocity  Splitting  to  the  Navier 
Stokes  Equations 

The  sound  speed  and  velocity  splitting  concept  can  be  applied  to  the  Navier 
Stokes  equations  because  of  the  special  form  of  the  viscous  Jacobi ans.  The 
Navier  Stokes  equations  in  strong  conservation  law  form  can  be  represented 


d%Q  +  dxE  -h  dyF  =  dxEvx  -|-  dxEvy  -f-  dyFvx  -j-  dyFvy 


(Al) 


where  Q,  E,  and  F  have  their  previous  definitions,  the  viscous  flux  terms  are 
given  by 


Evx  = 
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with 
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p  =  Pr-1(l-lT 


Beam  and  Warming  have  developed  a  general  class  of  integration  techniques 
for  equation  (A.1)  which  is  illustrated  here  for  Euler  implicit  time  differencing 


AQ"  -  At[«,AE"  +  SyAF’  -  },AE;,  -  SyAf;,] 

=  -  At[*,E”  +  V”  -  h(K,  +  K,)  -  1,(K,  +  *?,)] 


(A.3) 


■where  AQn  =  Qn+1  —  Qn  and  5  is  a  midpoint  operator.  In  this  algorithm  the 
cross  derivative  flux  terms  Evy  and  Fvx  are  treated  explicitly.  As  discussed  by 
Beam  and  Warming,  this  circumvents  the  need  to  include  an  implicit  viscous 
Jacobain  for  these  cross  terms,  and  it  also  maintains  unconditional  stablility 
for  the  scalar  model  equation. 

The  spacial  flux  terms  have  been  locally  linearized  in  time.  In  locally 
linearizing  these  terms  Beam  and  Warming  expand  the  viscous  terms  in  a 
Taylor  series  using  both  the  function  and  its  derivative,  as  for  example 
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Alternately,  these  terms  can  be  expanded  in  terms  of  Q  alone,  for  example 

dEvx 
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vx 
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where  is  now  a  differential  operator.  This  latter  from  will  be  used  here. 

Typically  p  is  not  linearized  in  time  and  the  Prandtl  number  is  treated  as  if 

it  where  a  constant  value  as  far  as  time  linerization  is  concerned.  For  a  term 
such  as  [idyll,  [i  Is  frozen  in  time  and  podvu  expanded  such  as 

fi0dy(pu/p)  =  fiodyftpotio/po)  —  {pou0/po2){p  —  Po )  +  Po~l[pu  —  pou0)] 
or  in  general 
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and 
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With  introduction  of  the  sound  speed  and  velocity  matrix  substitutions  for 
the  convection  Jacobian  matrices,  we  can  obtain  operators  Lx  and  Ly  that 
are  as  reducible  as  before  because  of  the  form  of  Avx  and  Bvy.  To  illustrate 
this  we  will  drop  the  numerical  dissipation  terms  as  before.  Then  Eq.(A.3) 
becomes 

L,L,AQ”  =  -  A ([«,£"  +  SyF"  -  },(K*  +  KJ  ~  lr(K*  +  ^)] 

(A.8) 
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Here  I  is  the  identity  matrix  and  cc  etc  are  elements  of  Ac  etc.  Keep  in 
mind  that  the  viscous  terms  have  been  second  order  accurate  differenced  to 
maintain  tridiagonal  structure  using  midpoint  operators  6. 

We  can  now  examine  the  structure  of  Lx  for  example,  and  show  that  the 
equations  uncouple  to  two  scalar  and  one  2x2  block  tridiagonal.  The  first 
row  of  Lx  is  clearly  uncoupled  and  can  be  solved  for  as  a  scalar  tridiagonal. 
Its  solution  is  then  used  to  update  the  right  hand  side  (RHS)  of  rows  2, 3,  and 
4.  With  the  first  column  of  the  matrices  brought  to  the  RHS,  the  third  row 
is  seen  to  be  uncoupled  and  can  thus  be  solved  as  another  scalar  tridiagonal. 
The  RHS  terms  of  the  second  and  fourth  rows  can  now  be  updated  and  what 


remains  is  a  2  x  2  block  tridiagonal  to  be  solved.  The  Ly  operator  inverts  in 
a  similar  fashion. 

Thus,  the  addition  of  the  viscous  terms  in  Cartesian  coordinates  does 
not  prevent  the  decoupling  technique  that  was  successfully  used  for  the 
inviscid  gasdynamic  equations.  Preliminary  work  shows  that,  using  the  local 
transformation  of  Section  HI,  the  same  decoupling  may  be  possible  for  the 
Navier  Stokes  equations  in  generalized  coordinates.  This  matter  is  being 
further  investigated. 


Appendix  B. Stability  of  Sound  Speed  and  Velocity  Splitting 

In  order  to  study  the  stablility  of  the  algorithm  proposed  in  Section  II, 
the  one  dimensional  Euler  equations  will  be  tested  for  linear  stability  using  a 
Fourier  (  or  matrix)  analysis. 

The  one  dimensional  Euler  equations  in  Cartesian  coordinates  are  given  by 

dtQ  +  dxE  =  0  (B.l) 


where 


Q  = 


'/>■ 

■  pu  - 

pu 

,  E  = 
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.  e . 

u(c  +  p). 

(B.2) 


A  frozen  coefficient  form  of  Eq.(B.l)  is  given  by 

dtQ  +  A*dxE  =  0 


(B.3) 


where  A  is  the  true  Jacobian,  [dE/dQ],  and  *  denotes  that  A  is  evaluated 
using  a  frozen  value  of  Q,  Q* .  Implicit  first  order  accurate  differencing  of 
Eq.(B.3)  is  given  in  delta  form  as 


I  +  hA*6^{Qn+l  -  Qn)  =  —  hA*6xQn 


(BA) 


In  our  sound  speed  -  velocity  splitting  scheme  we  replace  the  A*  matrix  of 
the  left  hand  side  operator  with  the  similarity  matrices  u*I  +  A*  to  obtain 


1  +  M.v! +Ae,)dw°+1  -  Q")  =  -hA'6,Q" 


(B.5) 


The  right  hand  side  matrix  A *  is  not  altered  in  our  procedure,  and  must  not 
be  changed  in  Eq.(B.5)  as  it  represents  the  correct  linearization  of  E  about 

Q *,  The  question  to  be  addressed  is  Eq.(B.5)  as  stable  as  Eq.(B.4)  now  that 
the  implicit  operator  as  been  altered? 

If  we  were  to  perform  a  stability  analysis  of  Eq.(B.4)  we  would  diagonalize 

A*  using  its  eigenvectors.  We  can  then  readily  perform  the  stability  analysis 
for  difference  equations  corresponding  to  scalar  partial  differential  equations. 
For  Eq.(B.5)  we  can  use  the  eigenvectors  of  A*  to  diagonalize  u  / -\ -Ae)  but  the 
right  hand  side  matrix  will  not  be  brought  into  diagonal  form.  Nevertheless, 
we  find  that  the  resulting  equations  can  still  be  analyzed. 


The  matrix  Ac  for  one  dimensional  flow  is  given  by 
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where 
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Multiplying  Eq.(B.5)  by  the  frozen-coefficient  inverse  of  the  eigenvectors 
(X*)— 1  gives 


[7  +  *(«’/ + aXV.](xT1W+1  -  Q")  =  -M^*r1A*x,«s(x*r1gn 

(B.10) 

or  in  nondelta  form  with  W  =  (X*)~lQ 
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However,  (X*)  1A*X*  —  u*/  —  A^  is  a  very  simple  matrix 
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Assuming  periodic  boundary  conditions  so  as  to  use  a  Fourier  stability 
analysis  (or  the  matrix  stablity  method  using  circulant  matrices),  the  finite 
difference  operator  6X  can  be  transformed  to 
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and  W  denotes  Fourier  transformation  of  W  (  i.e.  transformation  via  the 
eigenvectors  of  a  circulant  matrix). 

Applying  the  Fourier  method  to  Eq.(B.ll)  and  taking  the  inverse  of  the 
left  hand  side  operator  gives 
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or 
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The  amplification  matrix  associated  with  Eq.(B.14)  is  lower  triangular. 
Therefore,  its  eigenvalues  are  the  diagonal  elements,  and  these  clearly  have 
modulus  less  than  one  for  any  positive  value  of  h  =  At.  In  fact,  these  eigen¬ 
values  are  identical  to  those  obtained  with  the  standard  scheme  given  by 
Eq.(B.4).  Thus  the  necessary  condition  for  stability  is  met  for  any  values  of 
h.  A  sufficient  condition,  for  stability  requires  bounding  the  norm  of  the  3x3 
amplification  matrix.  Observing  the  3, 1  and  3, 2  elements  it  is  clear  that  the 
too  norm  can  be  unbounded  for  u  -*■  0  and  h  -+  oo,  however,  even  in  this 
norm  powers  of  the  amplification  matrix  will  be  quickly  bounded  unless  u  is 
identically  zero. 

It  is  interesting  to  note  that  if  one  were  to  form  Ae  using  the  Jacobian  of 
Ee,  the  resulting  algorithm  is  unconditionally  unstable.  In  this  case  Ac  is 
given  as 
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where 
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ac4  2  =  -^2  +  c2/7(7-  1) 

Proceding  in  a  similar  fashion  as  abo-ve,  the  resultant  representation  in  Fourier 
space  for  (B.5)  is 


The  maximum  eigenvalue  of  the  amplification  matrix  in  (B.16)  can  be  shown 
to  be  greater  than  one  for  all  At’s  greater  that  zero  and  the  necessary 
condition  for  stability  can  not  be  met.  In  actual  numerical  experimentation 
where  numerical  dissipation  is  added,  small  values  of  At  could  be  used  to 
obtain  stable  results. 


FIGURES 


Figure  1.  Arithmetic  operatiou  counts  of  block  tridiagonals  using 
L-U  decomposition. 

Figure  2a.  Overview  of  157  x  33  grid  generated  about  NACA  C~i'2 
airfoil  using  hyperbolic  grid  generator. 

Figure  2b.  Detail  near  body  of  157  x  33  grid. 

Figure  2c.  Pressure  distribution  on  NACA  0012  airfoil  (Ma 0  =  .80 
and  a  =  0.0) 

Figure  2d.  Residual  history  comparison  between  reduced  matrix 
form  and  standard  algorithms. 

Figure  3a.  Overview  of  249  x  50  grid  generated  about  NACA  0012 
airfoil  using  hyperbolic  grid  generator. 

Figure  3b.  Detail  near  body  of  249  x  50  grid. 

Figure  3c.  Pressure  distribution  on  NACA  0012  airfoil  (¥«,  =  .80 
and  a  =  0.0). 

Figure  3d.  Pressure  contours  near  airfoil. 

Figure  3e.  Residual  history  for  reduced  matrix  algorithm  with  scaled 
At. 

Figure  4a.  Overview  of  193  X  33  grid  generated  about  NACA  0012 
airfoil  using  hyperbolic  grid  generator. 

Figure  4b.  Detail  near  body  of  193  x  33  grid. 

Figure  4c.  Pressure  distribution  on  airfoil  {Moo  =  .80  and  a  = 
1.25). 

Figure  4d.  Pressure  contours  near  airfoil. 

Figure  4e.  Residual  history  for  reduced  matrix  algorithm  with  scaled 
At  and  pentadiagonal  inversions. 
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